# copied from Final_Code_Replication_Final/02_Analysis/00_Main_Text/Table_01/

# Code for Table 1: Incapacitation Mechanisms

rm(list = ls())

## ---------------------------------------
## Load Packages 
## ---------------------------------------
require('AER')
require('ivpack')
require('dplyr')
require('data.table')

## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))

## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)

## ---------------------------------------
## Mediation Analysis
## ---------------------------------------

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2)"
case.controls <-   "as.factor(any_drug_2) +  as.factor(any_violent_2) + as.factor(fire_arms_2) +  as.factor(any_rob_2) + as.factor(any_dui_2) + as.factor(prior_offender_2)"
demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + regis_before"

outc.1 <- "vote2012"
endo.1 <- "pti"
inst.1 <- "judgeiv"
med.1 <- "incar_f"

form.1m <- formula(paste(med.1, "~", endo.1, "+", time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
form.2m <- formula(paste(outc.1, "~", med.1, "+", endo.1, "+", time.controls, "+", demo.controls, "+", case.controls, "|", endo.1, "+", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))


## ---------------------------------------
# Row 1 - Nov. 2008 - Nov. 2012
## ---------------------------------------
m1a3m1 <- ivreg(form.1m, data = data0812)
m1a3m2 <- ivreg(form.2m, data = data0812)
m1a3m3 <- ivreg(form.3, data = data0812)

coef(m1a3m2)["pti"] + coef(m1a3m1)["pti"] * coef(m1a3m2)["incar_f"]
direct.effect <- coef(m1a3m2)["pti"]
total.effect <- coef(m1a3m3)["pti"]
indirect.effect <- total.effect - direct.effect
ie.se <- sqrt(vcov(m1a3m2)["pti", "pti"] + vcov(m1a3m3)["pti", "pti"])
te.se <-  sqrt(vcov(m1a3m3)["pti", "pti"])
de.se <-  sqrt(vcov(m1a3m2)["pti", "pti"])

main <- as.data.frame(cbind(direct.effect, de.se, indirect.effect, ie.se, total.effect, te.se))
names(main) <- c("Direct Effect (DE)", "DE Std Error", "Indirect Effect (IE)", "IE Std. Error", "Total Effect (TE)", "TE Std Error")
print("Table 1, Row 1. Nov. 2008 - Nov. 2012")
print(round(main, 3))

print('% Incapacitated...')
print(round((mean(data0812$incar_f, na.rm=T)*100),1))
print("% Pretrial incarcerated...")
print(round((mean(data0812$pti, na.rm=T)*100),1))
print("N...")
print(m1a3m2$n)

## ---------------------------------------
## Row 2 - Nov. 2011 - Nov. 2012
## ---------------------------------------
data0812 <- data0812[bail_date >= "2011-11-06", ]


m1a3m1 <- ivreg(form.1m, data = data0812)
m1a3m2 <- ivreg(form.2m, data = data0812)
m1a3m3 <- ivreg(form.3, data = data0812)

coef(m1a3m2)["pti"] + coef(m1a3m1)["pti"] * coef(m1a3m2)["incar_f"]
direct.effect <- coef(m1a3m2)["pti"]
total.effect <- coef(m1a3m3)["pti"]
indirect.effect <- total.effect - direct.effect
ie.se <- sqrt(vcov(m1a3m2)["pti", "pti"] + vcov(m1a3m3)["pti", "pti"])
te.se <-  sqrt(vcov(m1a3m3)["pti", "pti"])
de.se <-  sqrt(vcov(m1a3m2)["pti", "pti"])

cat("\nTable 1: Incapacitation Mechanism\n")
cat("Printing Table 1, Row 2. Nov. 2011 - Nov. 2012...\n")
main <- as.data.frame(cbind(direct.effect, de.se, indirect.effect, ie.se, total.effect, te.se))
names(main) <- c("Direct Effect (DE)", "DE Std Error", "Indirect Effect (IE)", "IE Std. Error", "Total Effect (TE)", "TE Std Error")

print(round(main, 3))

print("% Incapacitated...")
print(round((mean(data0812$incar_f, na.rm=T)*100),1))
print("% Pretrial incarcerated...")
print(round((mean(data0812$pti, na.rm=T)*100),1))

# observations
print("N...")
print(m1a3m2$n)

